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We use the Detrended Cross-Correlation Analysis (DCCA) to investigate the influence of sun activ- 
ity represented by sunspot numbers on one of the climate indicators, specifically rivers, represented 
by river flow fluctuation for Daugava, Holston, Nolichucky and French Broad rivers. The Multi- 
fractal Detrended Cross- Correlation Analysis (MF-DXA) shows that there exist some crossovers in 
the cross-correlation fluctuation function versus time scale of the river flow and sunspot series. One 
of these crossovers corresponds to the well-known cycle of solar activity demonstrating a universal 
property of the mentioned rivers. The scaling exponent given by DCCA for original series at in- 
termediate time scale, (12 — 24) < s < 130 months, is A = 1.17 ± 0.04 which is almost similar for 
all underlying rivers at Icrconfidence interval showing the second universal behavior of river runoffs. 
To remove the sinusoidal trends embedded in data sets, we apply the Singular Value Decomposi- 
tion (SVD) method. Our results show that there exists a long-range cross-correlation between the 
sunspot numbers and the underlying streamflow records. The magnitude of the scaling exponent 
and the corresponding cross-correlation exponent are A € (0.76,0.85) and 7x £ (0.30,0.48), re- 
spectively. Different values for scaling and cross-correlation exponents may be related to local and 
external factors such as topography, drainage network morphology, human activity and so on. Mul- 
tifractal cross-correlation analysis demonstrates that all underlying fluctuations have almost weak 
multifractal nature which is also a universal property for data series. In addition the empirical 
relation between scaling exponent derived by DCCA and Detrended Fluctuation Analysis (DFA), 
A « (/isun + /iriver)/2 is Confirmed. 



I. INTRODUCTION 

Recently, due to the developments in the area of com- 
plex systems as well as data measurements and data anal- 
ysis, one can find many opportunities for examination 
and interpretation of climate change which exhibit ir- 
regular systems [1-8]. It is well shown that the climate 
system is enforced by the well-defined seasonal periodic- 
ity, however the existence of unpredictable perturbation 
and chaotic functioning lead to extreme climate events. 
Indeed the climate is a dynamical system affected by 
tremendous factors and variables, such as solar activity 
which is represented by Sunspot numbers in this research 
[4-10]. All factors that control the trajectory of such 
mentioned systems have enormously large phase space 
and evolve as non-stationary processes, consequently we 
should explore it with stochastic tools to achieve reliable 
results. Nowadays, it has been clarified that a remark- 
ably wide variety of natural systems can be character- 
ized by long-rangeSuch cross-correlations address scien- 
tists toward fractal geometry of the underlying dynami- 
cal systems and can hopefully help us to predict future 
events. Existence and determination of power-law cross- 
correlations would help to promote our understanding of 
the corresponding dynamics and their future evolutions 
[11-14]. Beside, many events which controls earths cli- 
mate, water runoff records assigned by rivers and sun 
activity play a crucial and survival roles for human life. 
The runoff water fluctuations are excellent climate cri- 



teria because they integrate evapotranspiration (output) 
and precipitations (input) over large areas. It is well 
accepted that the prediction of water runoff is funda- 
mental for different aspects of social and economical rea- 
sons, ranging from the prediction of floods and droughts 
to planning of agricultural conditions. As a result of 
the periodicity in precipitation, river flow has also strong 
seasonal periodicity [12]. It is worth to note that unlike 
other climate components, water runoff may be directly 
influenced by human activity, like agriculture, drainage 
network morphology and so on, consequently makes it 
hard to distinguish the artificial and natural effects on 
the river flow data. Finding some or at least a universal 
behavior for different streamflow fluctuations as well as 
quantifying the impact of sun activity on various tem- 
poral and spatial scales of water runoff fluctuations can 
improve the recent hydrological models [15]. 

To this end, the statistical and multifractal analysis of 
river flows as well as influence of sun activity due to the 
interior and exterior chemical and physical properties of 
sun should be an important issue in the geophysical and 
hydrological systems. 

The streamflow of rivers and sun activity have been 
studied from various point of views such as: the probabil- 
ity distribution [16,17], correlation and fractal behaviors 
[18-23], connection between volatility and nonlinearity 
of fluctuations [11,12,24,25], scale invariancc for distribu- 
tion function [26] . In addition, sun activity have been in- 
vestigated by some methods in chaos theory [27] and also 



multifractal analysis [28,29], wavelet analysis [30], cross- 
correlation functions between monthly mean sunspot ar- 
eas and sunspot numbers [31,32], the relation between 
sunspot numbers fluctuation and number of flares, their 
evolution step [33,34], principal components and neural 
network methods to predict sunspots [35,36], sunspot ar- 
eas time series and solar irradiance reconstructions [37], 
magnetic and dynamic properties of sunspots at the pho- 
tospheric level [38] and the hydraulic-geometric similarity 
of river [39-41]. 

More recently, Pablo J. D. Mauas et. al., have investi- 
gated the solar forcing on climate, using the quantifica- 
tion of cross-correlation between the yearly sunspot num- 
bers, irradiance reconstruction and streamflow of Parana 
river [9,10]. On the other hand, the mechanisms for so- 
lar influence on the earth's climate has been clarified in 
detailed from various point of views in [42]. Q. Zhang et. 
al., have investigated the universal behavior of stream- 
flow records of the Pearl river [15]. 

After innovation of Hurst to propose the self-similar 
processes and its criteria, namely "Hurst exponent" 
[18-23], long-range correlated fluctuation behavior has 
also been reported for vast category of sciences, specif- 
ically the geophysical records (for more discussion see 
[18-23,43-46]). In the last decade, the modification pre- 
scription which is required for a full characterization of 
many data sets such as the runoff records, the various 
moments of the so-called fluctuation functions, have been 
introduced [18-23]. The effect of non-stationarity on the 
detrended fluctuation analysis has been investigated in 
[47,48]. 

Here we take a new approach and rely on the state- 
of-the-art algorithm to investigate the contribution of 
sinusoidal trends embedded in the data set as well as 
non-stationarity properties of the underlying series. We 
implement robust methods to explore the multifractal na- 
ture of cross-correlation between two important climate 
variables, the monthly streamflow of some rivers and sun 
activity represented by sunspot numbers (see Figure (1)), 
by using the novel approach in the fractal analysis, De- 
trended Cross-Correlation Analysis (DCCA) and its mul- 
tifractal modification, the Multifractal detrended Cross- 
Correlation Analysis (MF-DXA) [49,50]. We restrict this 
article to use the sunspot numbers as the solar activity 
indicator, since there are many large and reliable data 
sets which can be considered as solar influence on the 
climate. Due to the presence of the sinusoidal trends in 
both sunspot numbers and the runoff river fluctuations 
and based on previous researches, one cannot expect to 
find a unique scaling behavior for fluctuation functions 
in all time scales (see section HI), consequently, we have 
been motivated to use the well-known method, namely 
Singular Value Decomposition (SVD) method to exclude 
dominant trends in data set (see section II for more de- 
tails). So after, clean data set will be used in the DCCA 
and MF-DXA methods. 
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FIG. 1. Upper panel corresponds to the monthly sunspot 
number data set. The secular trend, obtained with a low-pass 
Fourier filter is shown as a thick line in the upper panel. Lower 
panels indicate observed flux fluctuations of Daugava, French 
Broad, Nolichucky and Holston rivers, respectively. The inset 
plot shows river flow for small scales. 

This paper is organized as follows: in section II, we 
describe the methods which are used to determine the 
cross-correlation of two non-stationary time series, the 
Detrended Cross- Correlation Analysis (DCCA), and in- 
vestigate the corresponding multifractal properties by us- 
ing the Multifractal Detrended Cross-Correlation Anal- 
ysis (MF-DXA). Section II will be continued by intro- 
ducing a method to eliminate trends from the original 
data set, the Singular Value Decomposition (SVD), and 
describing data used in this paper. In section HI, the 
multifractal cross-correlation of the underlying data sets 
will be examined. Section IV, will be devoted to the re- 
sults and summary. 



II. ANALYSIS TECHNIQUES AND DATA 
DESCRIPTION 



A. DCCA and MF-DXA 



Time series measured in the nature are usually af- 
fected by non-stationarities such as trends and artificial 
noises which must be well distinguished from the intrin- 
sic fluctuations of the series. In many cases also, in- 
trinsic fluctuations behave as non-stationary processes. 
Consequently, common methods in data analysis will 
be encountered with spurious or at least unreliable re- 
sults. One of the most famous and well-known approach 
used in many studies is Multifractal Detrended Fluc- 
tuations Analysis (MF-DFA) [51,52]. This method has 
been applied to various areas, such as economical time 
series [53-55,7,56], river flow [13] and sunspot fluctua- 
tions [57,58], cosmic microwave background radiations 
[59], music [60,61], plasma fluctuations [62]. 

For many reasons, wc arc interested in studying the 
mutual influence of two series in the presence of non- 
stationarities. Obviously, traditional methods for this in- 
vestigation become inaccurate procedures. Recently Jun 
et. al. have proposed an approach for analyzing cor- 
relation properties of a series by decomposing the origi- 
nal signal into its positive and negative fluctuation com- 
ponents [63]. Based on the previous study, Podobnik 
et. al. have modified the mentioned correlation method 
and improved it to explore the cross-correlation be- 
tween two non-stationary fluctuations, named Detrended 
Cross-Correlation Analysis (DCCA) [49] and its general- 
ized, the Multifractal Detrended Cross-Correlation Anal- 
ysis (MF-DXA) which also examine higher orders de- 
trended covariance [50]. 

As mentioned before, trends in data set may influence 
the accuracy of results. For reliable detection of the cross- 
correlations, it is essential to distinguish trends from the 
intrinsic fluctuations in data. Generally, trends embed- 
ded in measurements are of two types: Polynomial and 
Sinusoidal trends. Although the MF-DFA and MF-DXA 
methods eliminate the polynomial trends, the sinusoidal 
trends remain [47,48]. There are several robust meth- 
ods to eliminate the sinusoidal one such as Fourier De- 
trended Fluctuations Analysis (F-DFA) [62,64], which is 
actually a high-pass filter, and Singular Value Decompo- 
sition (SVD) [65,66]. One of the most disadvantage of 
the F-DFA method is the reduction of the size of under- 
lying data set. To apply the MF-DXA (see the following 
subsection) , we have to synchronize two underlying data 
set which may not be done using the F-DFA method. 
On the other hand, the SVD method promises to remain 
the length as well as synchronization [65,66]. Here in or- 
der to eliminate the effect of sinusoidal trends, we apply 
the Singular Value Decomposition (SVD) [65,66]. After 
trends elimination, we use the MF-DXA to analyze the 
cleaned data sets. 



One of the newly methods in analyzing two non- 
stationary time series is Detrended Cross-Correlation 
Analysis (DCCA) [49,63]. This method is a generaliza- 
tion of the DFA method in which only one time series was 
analyzed. Recently a generalized version of the DCCA 
method which is so-called MF-DXA, has been introduced 
[50]. Just same as the MF-DFA method, MF-DXA con- 
sists of the 4 steps (see [49,50,67] for more details): 
(I): Computing the profiles of the underlying data series, 
Xk and yk, as 

i 

fc=i 

i 

Y{t)^Y.^yk-{y)] t = l,...,N (1) 
fc=i 

the subtraction of mean is not compulsory. Since we are 
going to compare two different time series, we construct 
data sets with zero mean and unit variance using initial 
ones. 

(II) : Dividing each of the profile into Ng = int(A^/s) 
non-overlapping segments of equal lengths s, and then 
computing the fiuctuation function for each segments. 
In order to take the whole scries into account when the 
size of the data sets is not a multiple of considered time 
scale, s, we do the same procedure from the opposite end, 
consequently one finds 2Ns segments. 

1 " 

F{s,m) = - J^ {Y[{m - l)s + i] ~ yrn{i)} 
i—i 

x{X[{m-l)s + i]-x,„{i)} (2) 

for m = 1, A^s and: 
1 " 

F(s, m) - y{Y[N - (m - 1)^ + i] - y„,{i)} 

x{X[N - {m-l)s + i\- Xm{i)} (3) 

for m = Ng -|- 1, 2iVs, where Xm(i) and ym{i) are a fit- 
ting polynomial in segment mth. Usually, a linear func- 
tion is selected for fitting the function. If there is no 
trend in the data, a zeroth-ordcr fitting function might 
be enough [68]. 

(III) : Averaging the local fluctuation function over all 
the part, given by: 

F,^^) = \jfT.iF^'^^^)'\'"\ (4) 

L * rn=l J 

Generally, q can take any real value, except zero. For 
g = 0, equation (4) becomes: 



^o(s) = exp[-^^lnF(s,m)) 

\ " m=l / 

For q ~ 2, the standard DCCA is retrieved. 



(5) 



(IV): The final step is calculating the slope of the log- 
log plot of Fq{s) versus s which directly determines the 
scaling exponent X{q), as: 



(6) 



If both underlying series are equal then A(g) is noth- 
ing else but so-called generalized Hurst exponent, h{q). 
In the absence of sinusoidal trends embedded in data 
sets, if one finds no scaling behavior for the fluctuation 
function in equation (6) or at least, there does not exist 
any unique exponent for all scaling ranges then there ex- 
ists either short-range cross-correlation or not at all any 
cross-correlation. For a series of size N, the minimum 
number of windows will be Ns = 2 corresponding to the 
maximum value of s = int(7V/2). In addition, it has 
been demonstrated that to find the most correct value of 
scaling exponent by using DFA and DCCA methods, we 
should set s < {N/2), namely TVs > 2 [52]. 

To determine the slope of curve in the log-log plot of 
fluctuation function versus scale (equation (6)), we use 
Bayesian statistics [69] . We introduce measurements and 
model parameters as {X} : {Fq{s)} and {9} : {A(g)}, 
respectively. Based on the Bayesian theorem, the con- 
ditional probability of the model parameters given data 
set (observation) is so-called posterior probability and is 
given by: 



P{\{q)\X) = 



C{X\\{q))P{\{q)) 
J C{X\X{q))dX{q) 



(7) 



where the first term in the nominator of the right hand 
side is Likelihood and the second term contains all ini- 
tial constraints concerning model parameters, so-called 
prior distribution. This term expresses the degree of 
belief about the model. In the absence of every prior 
constraints, the posterior function, P{X{q)\X) is propor- 
tional to the Likelihood function. If there is no corre- 
lation between various measurements, consequently ac- 
cording to the central limit theorem. Likelihood function 
is given by a product of Gaussian functions as follows: 



C{X\X{q)) ^ C^p 



(8) 



where: 



xHX{q)) 



^„ [j;b..(^)" j"Thc.(g;A(g))]^ 

^obs.(s) 



Here Fobs.(s) and J^ThG.(s;A) are fluctuation function 
computed directly from the data set by using DFA or 



DCCA and determined by equation (6), respectively. 
Also, CTobs. (s) is the mean standard deviation, associated 
to -Fobs, (s)- Apparently, this Likelihood function to be 
maximum when for a value of the scaling exponent, X{q), 
reaches to its global minimum. The value of error-bar 
at 1(7 confidence interval of X{q) is determined by the 
likelihood function based on the following condition: 



68.3% 



C{X\X{q))dX{q) 



(10) 



Finally we report the best value of scaling exponent 
at la confidence interval according to A^'^_ for each mo- 
mcnt, g's. 

B. Singular Value Decomposition (SVD) 

Determining trends and construction proper detrend- 
ing operations are important step toward robust analysis, 
specially in climatic data analysis. As given by Z. Wu 
and his collaborators [70], there is no unique definition 
of trend and any proper algorithm for extracting it from 
underlying stationary as well as non-stationary data sets. 
In another aspect, the trend in a real world data series, 
non-stationary one, is an intrinsic function imposed by 
the nature on data set. To identify the trend on a data 
set, we can investigate the series in whole domain or on 
some specific span of domains. For linear and stationary 
data sets choosing the length of data set as domain of 
trend may be suit but for a real world data set which is 
non-stationary and nonlinear, we need more precise def- 
inition of trend. As of the importance of investigating 
trends and probably removing them from series, one can 
point out to two following aims: 

I) In some cases, there exists one or more crossover (time) 
scales, Sx, m the log-log plot of Fq{s) versus s (equa- 
tion (6)), segregating regimes with different scaling ex- 
ponents. These patterns demonstrates that underlying 
fluctuation has different correlation behavior in various 
values of scales [47,48,65,66,71]. 

II) In many cases, crossovers are produced by the em- 
bedded sinusoidal trends, e.g. seasonal trends in the cli- 
mate time series. Subsequently, to flnd scaling exponent 
of the intrinsic fluctuations, we should remove sinusoidal 
trends by using most robust detrending method so after, 
produced clean data set is used as an input for DCCA 
and MF-DXA methods [47,48,62]. 

In order to remove trends corresponding to the low fre- 
quency periodic behavior, we use Singular Value Decom- 
position method [72] instead of transformation a recorded 
data to the Fourier space using the method proposed in 
[73] (see also [47,48,71]). Using the SVD method not only 
we can track the influence of sinusoidal trends on the re- 
sults, but also the synchronization of two underlying data 
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FIG. 2. Power spectrum of the original sunspot numbers 
data set (solid line) and that of for a cleaned one using SVD 
method (filled symbols). Dashed lines correspond to the scal- 
ing function as for the same signals without sinusoidal 
trends. 
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FIG. 3. The DFA results of the original sunspot number 
data set and filtered one using SVD. The solid line corre- 
sponds to power fit of fiuctuation. The slope of power fits has 
been calculated by likelihood analysis (see the text). 



sets will be guaranteed [65,66]. In addition, we can de- 
termine over which scale, noises or trends have dominant 
contribution also the value of so-called crossover in the 
fluctuation function, in terms of the scale is computed by 
using DCCA method. [64-66,74,75]. After removing the 
dominant periodic functions, such as sinusoidal trends, 
we obtain the fluctuation exponent by direct application 
of the MF-DXA. 

The SVD method consists of the following steps: 

(I): Consider a data set which is superimposed with 
periodic trends {xi}; i = 1, N. Embed Xi with param- 
eters (d, r) where d and r are the embedding dimension 
and the time delay, respectively. The embedded data can 
be represented as a matrix T given by: 



where U^xd and V 



{N-{d-l}T}x(N-{d-l)T) 



arc the left and 



/ Xi Xi+r 



Xl+N-{d-l)T-l 



\ 



^i-\-N-{d-l)T-l 



(11) 



\ Xd Xd+T ■■■ a^d+7V-(d-l)T-l / 



as one can see, 1 < i < d. For a periodic embedded 
trend, we set p as the number of dominant frequency 
components of the form e'-"^'=*-', k ~ l...p in their power 
spectrum. For a fixed value of sample size, N, the max- 
imum value of d in the so-called trajectory matrix, F, 
equates to d < N - {d - 1)t + 1 [66,74,76]. 

(II): Matrix T will be decomposed to the two left and 
right orthogonal matrices as follows: 



USV^ 



(12) 



right orthogonal matrices, respectively. The diagonal el- 
ements of Sdx{N-{d-i)T) are the desired singular values, 
also known as eigenvalues. The Singular Value Decom- 
position of r is related to the eigendecomposition of the 
symmetric matrices F^r and TT\ as F^Fvi = Afvi and 
FF^Ui = Afui. The nonzero eigenvalues of F^F are the 
same as that of FF^'; and determine the rank of F. The 
eigenvalues are ordered such that Ai > A^+i > 0. The di- 
agonal elements of S will be constructed by the ordered 
eigenvalues (the others will be set to zero). The rank of 
F is equal to the number of nonzero eigenvalues. The 
columns of U and V arc constructed by the Ui and Vi, 
namely the eigenvectors mentioned before, corresponding 
to the ordered eigenvalues. 

By applying the SVD to the matrix F (i.e., F = USV^) 
we will get U, S and V. We consider the number of 
frequency components in the periodic trend to be p. 
Set the dominant 2p + 1 eigenvalues in the matrix S 
to zero and hereafter named as S*. The filtered ma- 
trix T* dx{N-{d-i)T) determined by F* = US*V^^ with 
elements r*^ . This in turn is mapped back on to a one- 
dimensional or filtered data given by: 



— ^ ij 



(13) 



where 1 < i < and 1 < j < iV - (d - 1)t. The p 
dominant cigcn-valucs and associating cigendecomposcd 
vectors, represent trends subspace subsequently, the re- 
maining (d — p) eigenvalues and corresponding eigenvec- 
tors demonstrate intrinsic fiuctuations subspace. In or- 
der to determine the value of p for a typical series such 



as monthly sunspot data set, at first, wc compute the 
power spectrum of monthly sunspot numbers. As shown 
in Figure (2), there arc at least two dominant sinusoidal 
trends embedded in the monthly sunspot numbers. The 
first one corresponds to the well-known sun activity and 
second period belongs to the so-called Schwabe cycle in- 
terval. In order to eliminate mentioned trends, we set 
p = 2 in the SVD algorithm and compute the power 
spectrum of the cleaned data again. Our expectation is 
that there must be no deviation from scaling function as 
v-!^ with /3 = 2h{q = 2) - 1 = 1.24 ± 0.02 as one can see 
in the lower plot in Figure (2). It is worth to note that by 
increasing p from its optimum value in the SVD method, 
probably some intrinsic statistical properties of underly- 
ing data set will be destroyed. After which we use the 
detrended sunspot series as an input data for common 
DFA method (see Figure (3)). As shown in Figure (3), 
all of the crossover time scales which are produced due to 
the competition between sinusoidal trends embedded in 
sunspot series and intrinsic fluctuations, were diminished 
and intrinsic fluctuations will be retrieved. Our result is 
in agreement with the previous result regarding to scaling 
behavior of sunspot based on MF-DFA method accom- 
panying by Fourier-Detrended Analysis, reported in [57]. 
Figure (3) conflrms that, SVD method could remove si- 
nusoidal trends. 

In the log-log plot of fluctuation function versus time 
scale given by DFA method, also one can find three 
crossovers (see Figure (3)). To determine their value, 
we define error function as: 

A(s) = 7[Fobs.(s)-i^Li„car(s)]' (14) 

for each q, where i^obs. (s) and FLinear(s) are the fluctua- 
tion functions for the original data and the filtered data 
produced by the SVD method, respectively. In Figure (4) 
we plot A as a function of s for the sunspot numbers fluc- 
tuations. The first crossovers occurs at six ~ [12 — 18] 
months corresponding to the annual period. The second 
crossovers is equal to S2x ^ [130 — 170] months which is 
related to the well-known solar activity period. We can- 
not determine the value of third crossover with the ac- 
ceptable uncertainty by using DFA method due to small 
size of current sunspot series. 

C. Data Description 

We use the up-to-date monthly sunspot numbers {Sn) 
data series released by National Oceanic and Atmo- 
spheric Administration (NOAA) [77] and the Sunspot 
Index Data Center (SIDC) [78]. The monthly flow fluctu- 
ations of four famous rivers, namely Daugava at Latvia, 
Holston near Damascus, Nolichucky at Embreeville and 
French Broad at Asheville were collected from National 
Water Information System [79]. The original Daugava 
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FIG. 4. A function versus time scale. 

river data source is Latvian Environmental Geological 
and Meteorological Agency database [80] . The runoff di- 
mension of the underlying rivers is m'^/s (see Figure (1)). 
The river flow fluctuations have been measured indepen- 
dently in corresponding area. All mentioned rivers are 
mixed feeding from rain, snowmelt water and groundwa- 
ter. Table (I) reports some characteristics of river flow 
fluctuations used in this study. These data sets are al- 
most long in length and most available series. In Figure 
(5) we plot underlying streamflows and sunspot numbers 
detrended data sets. In this figure, we use SVD method 
explained in the previous subsection for detrending data 
sets, also the mean and variance of series equate to zero 
and unity, respectively. According to Figure (5), one 
may find a remarkable correlation for streamflows and 
sunspot numbers. To quantify this correlation we use 
the Pearson's correlation coefficients for X and Y de- 
fined by: r ~ 1'''"^^^^^^'^"^'^^^ . The value of r for de- 
^ ^ V(x-(x»(y-(y» 

trended river flows and sunspot numbers correspond to 

^Daugava = +0.44, THolston = +0.94, rNolichucky = +0.88 
and rprcnch-Broad = +0.94. 

III. MF-DXA OF SUNSPOTS AND RIVER FLOW 
FLUCTUATIONS 

To examine the multifractal properties and cross- 
correlation of sunspot numbers and river flow fluctua- 
tions, we apply the DCCA and MF-DXA methods. As 
mentioned in the previous section, we have to synchronize 
two time series of interests. The length of the river flow 
fluctuations may vary from 600 to 1300 months for the 
studied rivers, Daugava, Nolichucky, Holston and French 
Broad. On the contrary, there is more information avail- 
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FIG. 5. The detrended fluctuations of rivers (dashed 
line) and the corresponding sunspot (soUd line) time se- 
ries. Here we use a SVD method to trend series. The 
mean and variance of each series became zero and unity, 
respectively. The Pearson's correlation coefficients are 
^Daugava = +0.44, rHoiston = +0.94, rNoiichucky = +0.88 and 

''FronchBroad = +0.94. 

able for sunspot databases, ranging from daily to annu- 
ally data sets. To find reliable results given by DCCA, 
we synchronized the sunspot monthly series to all four 
river's data sets. 

Before applying any detrending program we apply 
DFA, DCCA and MF-DXA method. Log-log plot of F{s) 
versus s for original streamflow and sunspot data sets 
by using DFA method (Figure (6)) eonfirms that all un- 
derlying runoff rivers behave as the anti-persistent long- 
range correlated series for time scale s < 12 months, 
the value of Hurst exponent at this scale for Daugava, 
Holston, Nolichucky, French Broad and sunspot corre- 
spond to ^Daugava = h{2) - 1 = 0.21 ± 0.02, i?Holston = 

/i(2)-l = 0.22+0.02, T/Noiichucky = /i(2)-l = 0.21+0.02, 

-ffprench- Broad = /j(2) — 1 = 0.18 + 0.03 and i?Sunspot = 

h{2) -1^ 0.11 + 0.03, respectively. 

Figure (7) shows the fluctuation function given by mul- 
tifractal DCCA method for each river flow fluctuations 
versus sunspot numbers without applying any detrend- 
ing procedure. Some crossover time scales in fluctua- 
tion fimctions are detected in Figure (7). One of those 
crossovers is known as cycle of solar activity. In addition 
the scaling behavior of fluctuation function at interme- 
diate time scales, namely [12 — 24] < s < 130 months, is 
similar at Icr confidence interval for q = 2 for all rivers 
used in this paper (see Figure (8)). The scaling exponent 
at this regime is X{q = 2) = 1.17 + 0.04. As shown in 
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FIG. 6. Fluctuation function versus scale in log-log plot 
computed by DFA method for original data sets. To make 
more sense we shifted the vertical scale. 

TABLE I. The main characteristics of runoff rivers used in 
this paper. 



River 


Discharge 
(m^/sec) 


Series Length 


Drainage 
area(km^) 


Location 


French Broad 


2093 


1896.1 - 2005.12 


2448 


35°36' N 
82° 34' W 


Daugava 


601 


1953.1 - 2002.12 


87900 


56° 57' N 
24° 6' E 


Holston 


479 


1931.1 - 2005.12 


784.8 


36° 39' N 
81°50' W 


Nolichucky 


1379 


1921.1 - 2005.12 


2085 


36° 10' N 
82° 27' W 



Rcf. [48], trends become dominant at intermediate regi- 
men while at small and very large time scales, the intrin- 
sic fluctuations to be dominated, therefor, we conclude 
that at [12 — 24] < s < 130 months the cross-correlation 
of water runoff and corresponding sunspot numbers de- 
termined by DCCA method for (7 = 2, is characterized 
as the universal behavior. For s < [12 — 24] and s > 130 
months the slopes of fluctuation functions are different 
for various rivers and implies that at mentioned time 
scales the local effects such as morphology, human ac- 
tivities, various drainage areas become dominant [15,81]. 
Moreover, cycle of solar activity represented by sunspot 
is one of most robust effect which affects on streamflow 
fluctuations. 

Figure (9) shows the power spectrum of original and 
detrended mentioned rivers data sets and synchronized 
sunspot numbers series. The spectral density of French 
Broad and Nolichucky rivers behave in similar way 
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FIG. 7. Log-log plot of fluctuation functions given by 
the MF-DXA [FDCCAi^sy) method for various values of 
q for four rivers' flow mentioned in the manuscript ver- 
sus sunspot numbers as a function of time scale. Here 
Sraax = [size of series/2]. We shifted the vertical scale to 
make more sense. 
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FIG. 8. Log-log plot of fluctuation functions given by 
DCC A{F DC ca{s)) method for four streamflows fluctuations 
versus sunspot numbers as a function s at intermediate time 
scale, 12 — 24 < s < 130 months. For clearly, we shifted the 
vertical scale. Solid lines correspond to the power flt of fluc- 
tuation function and slopes have been computed by likelihood 
analysis. 



around the well-known solar activity. This outcome 
is also confirmed in Figure (7), since there is a sharp 
crossover in the fluctuation function versus s, in the MF- 
DXA method for French Broad and Nolichucky rivers. 
This may demonstrate that the contribution of solar ac- 
tivity depends on the geographical properties of river. 

Now, to explore the existence of cross-correlation in 
data set without sinusoidal trends from the DCCA and 
MF-DXA method we detrend each of the above synchro- 
nized data set by using the SVD method. Applying the 
DCCA method, we find out that there exists an strong 
long-range cross-correlation between sun activity repre- 
sented by sunspot numbers and river flow fluctuations. 
Figure (10) shows the log-log plot of fluctuation functions 
introduced by DFA {FupAis)) for a typical river flow, 
namely Daugava and sunspot fluctuations in the same 
time interval as well as the same function given by DCCA 
{Fdcca[s)) for those synchronized ones and q = 2. Fig- 
ure (11) indicates the results of DCCA method for three 
other rivers with their synchronized sunspot series. The 
scaling behavior of these functions represented in Figures 
(10) and (11) conflrms that there exists cross-correlation 
between such river flow fluctuations and sun activity. To 
explore the multifractal nature of underlying data sets, 
we apply MF-DFA [52,57] and MF-DXA [50] methods. 
Figure (12) shows generalized Hurst exponent and X{q) 
as a function of q for introduced series. Filled square 
and circle symbols show generalized Hurst exponent for 
river and corresponding synchronized sunspot data set, 



respectively. Filled triangle symbol indicates generalized 
scaling exponent, A(g). Weakly ^-dependency of scaling 
exponent implies the rivers and sunspot numbers have 
almost multifractal behavior. For all studied rivers, the 
segments with fluctuations in the MF-DXA method near 
the mean value are larger than that of far from the mean 
as shown in Figure (12), namely the value of X{q) for 
g < is almost larger than that of for g > 0. This 
shows the statistics of small fluctuations in water runoff 
arc bigger than large fluctuations, rare events. Accord- 
ing to Figure (12), the value of X{q) for g > remains 
almost constant, this indicates that the rare events in 
the cross-correlation function behave in similar way. On 
the other hand, the scaling behavior of fluctuation func- 
tion, Fq{s) in the MF-DXA method confirms the various 
type of fluctuations in different time scales have the same 
fractal features. This may be useful to extend the frac- 
tal properties of cross-correlation function which are de- 
rived in the short time scales to larger one in the absence 
of large time series. In addition, according to Figure 
(12), for g = 2, we obtain the empirical approximation, 
A(q = 2) « [/i,u„(<7 = 2) + /iHvor((Z = 2)]/2. 

According to auto-correlation function given by: 

C7(r) = {[xii + r) - {x)][x{i) ~ (x)]) - r"^ (15) 

we can introduce the cross-correlation function for so- 
called long-range cross-correlation behavior as: 

Cx (r) = {[x{t + r) - {x)M^) - {y)]) ^ r'^^ (16) 




FIG. 9. Spectral density of river flow fluctuations, de- (-Fdfa(s)) and DCCA {Fdcca{s)) methods for detrended 

trended data set and synchronizing sunspot numbers. We sunspot numbers and Daugava river flow as a function of time 

shifted the vertical scale for spectral density to make more scale. Where we choose q = 2. The slope of power fits have 

sense. been derived by likelihood statistics. 



where 7 and 7x are the auto-corrclation and cross- 
correlation exponents, respectively. Very often, direct 
calculation of these exponents are not recommended due 
to the non-stationarities and trends superimposed on the 
collected data. One of the reliable and proper statistical 
methods to calculate auto-correlation exponent is DFA 
method, namely 7 = 2- 2h{q = 2) [13,71]. Recently B. 
Podobniket. al. have demonstrated the relation between 
cross-correlation exponent, 7x and scaling exponent de- 
rived by equation (6) according to 7x = 2 — 2A(g = 2) 
[49] . The magnitude of cross-correlation and scaling ex- 
ponent derived by DFA and DCCA arc reported in Table 
II. 

TABLE II. Values of the scaling and cross-correlation ex- 
ponents of detrended sunspot numbers and river flow fluc- 
tuations using the MF-DXA method for four rivers as well 
as scaling exponent given by the DFA method, h{q), of each 
synchronized data set for q = 2. 



River 



Nolichucky 



French Broad 



Holston 



Daugava 



^sun 



spot 



0.93 ±0.01 



1.11 ±0.01 



1.05 ±0.01 



1.14 ±0.01 



0.70 ±0.01 



0.63 ±0.01 



0.60 ±0.01 



0.61 ±0.01 



A 



0.85 ±0.01 



0.84 ±0.01 



0.78 ±0.01 



0.76 ±0.01 



7x 



0.30 ± 0.02 



0.32 ±0.02 



0.44 ± 0.02 



0.48 ± 0.02 



To quantify the impact of El Nino (ENSO) phe- 
nomenon, we used the El Niiio 3 index reported in [82] 
since 1950. Applying SVD method, the sinusoidal trends 
have been diminished, thereafter, cleaned data sets used 



as inputs for DCCA algorithm. Figure (13) indicates the 
results computed by DCCA method for the detrended 
data sets. There is a crossover time scale in the fluctu- 
ation function as a function of scale for all underlying 
rives. The value of crossover equates to Sx ~ 60 months 
which is so-called the period of El Niiio phenomenon. 
The value of scaling exponents and cross-correlation ex- 
ponents reported in Table III, confirm that on s < 60 
months, there exists a cross-correlation between ENSO 
phenomenon and rivers. For time scale larger than s > 60 
months, this cross-correlation becomes ignorable [9,10]. 
We also compare the impact of ENSO phenomenon and 
sun activity on the streamflow of rives. To this end, we 
should synchronize the El Nifio, sunspot and river flow 
fluctuations. The value of scaling exponents and cross- 
correlation exponents for sunspot-river have also been 
reported in Table III. We find that the contribution of 
sun activity represented by sunspot is almost larger than 
ENSO phenomenon on the mentioned runoff rivers dur- 
ing mentioned period. 



IV. DISCUSSION AND CONCLUSION 

There are many motivations such as complexity of fluc- 
tuations in our environments which lead to examine them 
using robust and novel methods from complex systems 
and statistical physics point of views. Knowledge of nat- 
ural variabilities are necessary to manage the energy re- 
sources and prevent disasters from social and econom- 
ical point of views. In this paper, we analyzed three 
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FIG. 11. Fluctuation function given by DCCA method for 
detrended Holston, French Broad and Nolichucky rivers with 
their synchronized sunspot series in log-log scale versus time 
for q = 2. The maximum value of s for each series equates to 
the integer part of half size of the corresponding river flow. 

TABLE III. Values of the scaling and cross-correlation ex- 
ponents of synchronized El Nifio 3 index and river flow fluc- 
tuations as well as for sunspot and river in the same period, 
using the DCCA method. 



River 



Nolichucky 



French Broad 



Holston 



Daugava 



A 



ENSO 



0.89 ±0.03 



0.89 ±0.03 



0.85 ±0.03 



0.74 ± 0.03 



0.22 ±0.06 



0.22 ±0.06 



0.30 ± 0.06 



0.52 ±0.06 



nspot 



0.94 ±0.01 



0.98 ±0.02 



0.90 ±0.01 



0.77 ±0.01 



7 



sunspot 



0.12 ±0.02 



0.04 ± 0.02 



0.20 ±0.02 



0.46 ± 0.02 



important fluctuations in the nature, namely sun activ- 
ity illustrated by sunspot numbers, El Niiio phenomenon 
and streamflovir of rivers. However, many climate indica- 
tors such as river flow fluctuations are afl'ected by dom- 
inant seasonal trends, but recent researches have con- 
flrmed that, there are many variables causing the evo- 
lution of environmental conditions behave as a complex 
systems. Subsequently, applying the common methods in 
data analysis give incorrect or at least unreliable results. 
On the other hand to infer valuable results, the following 
necessary conditions should be satisfled: 

i) The length of underlying time series should be large 
enough. 

ii) The contribution of superimposed trends and noises 
on the recorded data must be small enough or at least 
distinguishable . 

Unfortunately above necessary conditions cannot be met 
in some practical measurements. To solve these prob- 
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FIG. 12. Generalized Hurst exponent, h{q) computed by 
MF-DFA (filled circle and square symbols) and A(q) derived 
by equation (6) in MF-DXA method (filled triangle symbol). 
Filled circle and square correspond to sunspot numbers and 
mentioned runoff river fluctuations, respectively. 



lems, we rely on the robust methods in data analysis 
to explore the mutual effect of sunspot and runoff wa- 
ter fluctuations. Based on the mentioned motivation, we 
apply the most recent method, Multifractal Detrended 
Cross-Correlation Analysis (MF-DXA), to examine the 
cross-correlation and fractal properties of sunspot num- 
bers and river flow fluctuations for some most famous 
and available rivers. 

Based on previous researches, as we expect, there are 
many sinusoidal trends embedded in the underlying sig- 
nals. Unfortunately, these trends cannot be removed 
by common detrended procedures in DFA and DCCA 
analysis [47,48]. These trends may cause some spurious 
crossovers in the fluctuations function versus time scale 
related to the competition between noise and trends. Ac- 
cording to Figure (4), there are some crossovers in the 
fluctuation function of the sunspot numbers versus time 
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FIG. 13. Fluctuation function versus scale in log-log plot 
computed by DCCA method for detrended river streamflows 
and El Nino index. To make more sense we shifted the vertical 
scale. 

scale, s. The first one corresponds to the annually pe- 
riod, the second is equal to the so-called solar activity, 
11 years, while the last one may indicate the well-known 
Gleissberg period, but we must point out that this value 
with small uncertainty cannot be determined by this cur- 
rent data, because the size of data set is small. Here to 
eliminate those trends, we use the Singular Value Decom- 
position (SVD) method. 

Results given by DFA method for original runoff rivers 
and sunspot numbers, demonstrated that all underlying 
runoff river behave as the anti-persistent long-range cor- 
related series for time scale s < 12 months at la confi- 
dence interval (Figure (6)). On the other hand, the value 
of Hurst exponents for series without sinusoidal trends 
exhibit runoff rivers and sunspot numbers are long-range 
correlated process (see Table H). 

MF-DXA analysis for non-detrended underlying data 
sets implied universal scaling exponent of fluctuation 
function for g = 2 for all underlying rivers flow at in- 
termediate time scales, [12 — 24] < s < 130 months, 
namely A = 1.17 ± 0.04 (see Figure (8)). This find- 
ing has recently been reported in Ref. [81] from power 
spectrum analysis. On the other hand there is a simi- 
lar crossover in cross-correlation fluctuation functions at 
Sx ^ 130 months which is known as cycle of solar activ- 
ity, for all studied rivers versus sunspot [13]. 

Data without sinusoidal trends were used as the inputs 
for DCCA and its generalized MF-DXA methods. Due 
to the scaling behavior of fluctuations function (equation 
(6)) versus time scale, s, we concluded that there exists a 
significant cross-correlation between sun activity and wa- 



ter runoff river fluctuations. This cross-correlation con- 
firms that the influence of other reasons on the river flow 
fluctuations such as human activity, drainage network 
morphology, land use patterns and topography may be 
ignored or at least cannot be distinguished form major 
effect which here is sun activity [81]. It must point out 
that, we know that the amount of radiation given off 
by the Sun (solar irradiance) is greatest when there are 
lots of sunspots. As discussed in detail in Ref. [9], the 
higher value of galactic cosmic ray, the higher the cloud 
cover which increases water resource of river, so we ex- 
pect streamflow to be almost cross-correlated with so- 
lar activity. Other explanation are also investigated in 
[83,84]. 

The cross-correlation exponent, 7x , of underlying data 
set in the presence of non-stationarities and trends have 
been determined by DCCA method. According to the 
relation between scaling exponent, A(g = 2), and cross- 
correlation exponent, 7x, namely 7x = 2 — 2X{q = 2), 
one flnd that for long-range cross-correlated signals, the 
value of 7x goes to small value demonstrating the cross- 
correlation function decreases more slowly and statisti- 
cally two underlying data sets tend their present situa- 
tion in their cross-correlation behavior. According to the 
values reported in Table II and Figures (10) and (11), 
there exists a long-range cross-correlation between rivers 
flow fluctuations used in this paper and sun activity in- 
dicated by sunspot numbers. 

Weak g-dependency of the generalized Hurst, h{q), and 
scaling, A(g), exponents, demonstrates that the origi- 
nal underlying series and also cross-correlation between 
sunspot numbers and each river behave as almost multi- 
fractal processes, demonstrating another universal char- 
acteristics. The value of A(g) for g > remains almost 
constant, this indicates that the rare events in the cross- 
correlation function behave in similar way. The simi- 
larity in behavior of fluctuation function, Fq{s) in the 
MF-DXA method demonstrates the various type of fluc- 
tuations in different time scales have the same fractal 
features. Based on this result, this may be useful to ex- 
tend the fractal properties of cross-correlation function 
derived by using short time scales to larger one without 
relying on the large time series. In addition, the rare 
events in streamflow and sunspot cross-correlation anal- 
ysis, are not affected by local characteristics. 

According to our results, the empirical relation be- 
tween X(q = 2) and Hurst exponent, X{q = 2) w 
[hsun{q = 2) + Kivci{q = 2)]/2, has also been confirmed. 

Since the value of cross-correlation of Daugava river 
and sunspot numbers is less than other rivers, one can 
conclude that, beside the sun as a main resource of en- 
ergy in the nature, the river's geographical situation and 
the source of rainfall for that river may have reasonable 
impact on runoff water fiuctuation. Recently, even an 
opposite behavior between solar activity and rainfall fluc- 
tuation in equatorial east Africa has been noticed in the 



recent report of the East Africa [85] explained in Ref. 
[9]. In addition, the value of A for French Broad and 
Nolichucky rivers are larger than two remained runoff 
rivers. Figure (9) also indicated, the power spectrum 
of French Broad and Nolichucky rivers behave in simi- 
lar way around the well-known solar activity period. We 
also found a sharp crossover in the fluctuation functions 
derived by MF-DXA method for mentioned rivers in Fig- 
ure (7). Since, the geographical positions of French Broad 
and Nolichucky rivers are close together, subsequently as 
expressed before, the impact of sun activity represented 
by sunspot, depends on the geographical properties of 
rivers. 

One of the most important phenomenon which can af- 
fect the runoff rivers is El Nifio (ENSO). To compare 
the contribution of ENSO phenomenon and sun activity, 
we used El Nifio index 3. The results derived by DCCA 
have been indicated in Figure (13) and also reported in 
Table III. There exists a crossover time scale in the log- 
log plot of fluctuation function versus scale. The value 
of this crossover time scale is Sx ^60 months regard- 
ing to well-known ENSO period. ENSO phenomenon is 
cross-correlated by streamflow fluctuations for s < 60 
months, while for larger than this period, one can ignore 
this cross-correlation behavior. By comparison, the ef- 
fect of ENSO phenomenon and sun activity on the river 
flows, we found that the contribution of sunspot is almost 
larger than ENSO index since 1950. 

There are many methods to predict the sun activity 
[35,36], so based on our current results which is demon- 
strated the cross-correlation between detrended sunspots 
and water runoff, one may use those results to predict 
river flow fluctuations. 

Finally, one must point out that it should be inter- 
esting to extend the present analysis to a various set of 
runoff water fluctuations and other sun activity indica- 
tors such as solar irradiance, galactic cosmic rays flux 
and geomagnetic activity to find whether the values for 
the main parameters used in this analysis have a more 
universal validity of the solar influence on the streamflow 
fluctuations assigned as a climate variable in this paper. 
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